# This file esimtates the difference in proportions at all bandwidths from 1 to 183

load("mock_datas_nolag.rdata")
load("mock_ns_nolag.rdata")

## commented out for replication files 

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("reshape2")
#  install.packages("dplyr")

library(xtable)
library(ggplot2)
library(rmeta)
library(reshape2)
library(dplyr)

# write function for two samples proportion test

pr_mean_se <- function(var, by){
  p_1 <- var[by == 1]
  p_0 <- var[by == 0]
  mean_1 <- summary(p_1)[4] 
  mean_0 <- summary(p_0)[4]
  pr_mean <- mean_1 - mean_0
  pr_se   <- sqrt( mean_1*(1-mean_1) / length(p_1) + 
                     mean_0*(1-mean_0) / length(p_0) )
  p       <- 1-pnorm(pr_mean/pr_se) 
  n_under <- length(p_0)
  n_over  <- length(p_1)
  return(c(pr_mean, p, pr_se, mean_0, n_under, n_over))
}

# create data frame for diff-in-prop estimates in close windows 
close_ests <- data_frame(close_ests = rep(NA, 183),
                         ps         = rep(NA, 183),
                         close_se   = rep(NA, 183),
                         close_tu   = rep(NA, 183),
                         n_under    = rep(NA, 183),
                         n_over     = rep(NA, 183))

data <- left_join(data_list[[49]],
                  data_n_list[[49]],
                  by = c("days", "year"))

for(i in 1:183){
  data_close <- 
    data %>%
    filter(abs(days) < i) %>%
    mutate(voted     = round(par_vote*n),
           abstained = round((1-par_vote)*n),
           treated   = days > 0) %>% 
    select(days, treated, voted, abstained)
  
  data_voted <- 
    as.data.frame(lapply(data_close, 
                         function(x,p) rep(x,p), 
                         data_close[["voted"]])) %>%
    mutate(voted = 1)
  
  data_abstained <- 
    as.data.frame(lapply(data_close, 
                         function(x,p) rep(x,p), 
                         data_close[["abstained"]])) %>%
    mutate(voted = 0)
  
  data_close <- 
    rbind(data_voted, data_abstained)
  
  test <- with(data_close, pr_mean_se(voted, treated))
  
  close_ests[i, 1] <- test[1]
  close_ests[i, 2] <- test[2]
  close_ests[i, 3] <- test[3]
  close_ests[i, 4] <- test[4]
  close_ests[i, 5] <- test[5]
  close_ests[i, 6] <- test[6]
  print(i)
}

close_ests$place <- 1:183

local_effect <- 
  ggplot(close_ests, 
         aes(x = place,
             y = 100*close_ests)) + 
  geom_line() + 
  geom_line(aes(x = place,
                y = 100*(close_ests + qnorm(0.005)*close_se)),
            alpha = 0.2, linetype = "dashed") + 
  geom_line(aes(x = place,
                y = 100*(close_ests + qnorm(0.995)*close_se)),
            alpha = 0.2, linetype = "dashed") +
  geom_line(aes(x = place,
                y = 100*(close_ests + qnorm(0.025)*close_se)),
            alpha = 0.6) + 
  geom_line(aes(x = place,
                y = 100*(close_ests + qnorm(0.975)*close_se)),
            alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  theme_classic() +
  xlab("Bandwidth in days") + 
  ylab("Effect on turnout")  +
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14))  

table_closeests <- 
  matrix(NA, nrow = 6, ncol = 5)

select <- c(1,2,3,7,30)
table_closeests[1,] <- round(100*close_ests$close_ests[select],2)
table_closeests[3,] <- round(100*close_ests$close_tu[select],2)
table_closeests[4,] <- round(close_ests$n_under[select],2)
table_closeests[5,] <- round(close_ests$n_over[select],2)

for (i in 1:5){
  table_closeests[2,i] <- paste("[", 
                                round(100*(close_ests$close_ests[select[i]]+qnorm(0.025)*close_ests$close_se[select[i]]),2),
                                "; ",
                                round(100*(close_ests$close_ests[select[i]]+qnorm(0.975)*close_ests$close_se[select[i]]),2),
                                "]", sep ="")

  prop_above <- with(close_ests, (n_over[select[i]] / (n_over[select[i]] + n_under[select[i]])))
  
  table_closeests[6, i] <- paste("[", 
                                 100 * round(prop_above + qnorm(0.025) * 
                                   sqrt(((1 - prop_above) * prop_above) / 
                                          (close_ests$n_over[select[i]] + close_ests$n_under[select[i]])), 4),
                                 "; ",
                                 100 * round(prop_above + qnorm(0.975) * 
                                   sqrt(((1 - prop_above) * prop_above) / 
                                          (close_ests$n_over[select[i]] + close_ests$n_under[select[i]])), 4),
                                 "]", sep ="")
}

colnames(table_closeests) <- c("1 day", "2 days", "3 days", "1 week", "1 month")
rownames(table_closeests) <- c("Effect on turnout", "95% CI", "Turnout below cutoff", 
                               "N below cutoff", "N above cutoff", "CI for distribution of N")

xtable(table_closeests)
